
/* --------------------------------------------------------------------
   To compile:
   cc -g -c polylib.c -I ./
   cc -g -c polylib_test.c -I ./
   cc -g -o polytest polylib_test.o polylib.o -lm 
 * --------------------------------------------------------------------*/

/* -------------------------------------------------------------------
   This is a routine to test the integration, differentiation and
   interpolation routines in the polylib.c. 

   First, it performs the integral

      /1      alpha   beta  alpha,beta
     |   (1-x)   (1+x)     P (x)       dx  = 0
     /-1                    n

   for all   -0.5 <= alpha <= 5   (increments of 0.5)
             -0.5 <= beta  <= 5   (increments of 0.5)

   using np points where
	      NPLOWER <= np <= NPUPPER
                2     <= n  <= 2*np - delta

   delta = 1 (gauss), 2(radau), 3(lobatto).
   The integral is evaluated and if it is larger that EPS then the 
   value of alpha,beta,np,n and the integral is printed to the screen.

   After every alpha value the statement
       "finished checking all beta values for alpha = #"
   is printed

   The routine then evaluates the derivate of 

          d   n      n-1
	  -- x  = n x 
	  dx

   for all   -0.5 <= alpha <= 5   (increments of 0.5)
             -0.5 <= beta  <= 5   (increments of 0.5)

   using np points where
	      NPLOWER <= np <= NPUPPER
                2     <= n  <= np - 1

   The error is check in a pointwise sense and if it is larger than 
   EPS then the value of alpha,beta,np,n and the error is printed to 
   the screen. After every alpha value the statement
       "finished checking all beta values for alpha = #"
   is printed

   Finally the routine  evaluates the interpolation of

             n      n
	    z  to  x

   where z are the quadrature zeros and x are the equispaced points
       
                  2*i
        x    =   -----   - 1.0    (0 <= i <= np-1)
	 i       (np-1)
            
	  
   for all   -0.5 <= alpha <= 5   (increments of 0.5)
             -0.5 <= beta  <= 5   (increments of 0.5)

   using np points where
	      NPLOWER <= np <= NPUPPER
                2     <= n  <= np - 1

   The error is check in a pointwise sense and if it is larger than 
   EPS then the value of alpha,beta,np,n and the error is printed to 
   the screen. After every alpha value the statement
      "finished checking all beta values for alpha = #"
   is printed

   The above checks are performed for all the Gauss, Gauss-Radau and
   Gauss-Lobatto points. If you want to disable any routine then set
      GAUSS_INT, GAUSS_RADAU_INT, GAUSS_LOBATTO_INT = 0 
   for the integration rouintes
      GAUSS_DIFF,GAUSS_RADAU_DIFF, GAUSS_LOBATTO_DIFF = 0
   for the differentiation routines 
      GAUSS_INTERP,GAUSS_RADAU_INTERP, GAUSS_LOBATTO_INTERP = 0
   for the interpolation routines.
   ------------------------------------------------------------------*/
   
#define NPLOWER  5
#define NPUPPER 15
#define EPS  5e-12

#define GAUSS_INT            1
#define GAUSS_RADAUM_INT     1
#define GAUSS_RADAUP_INT     1
#define GAUSS_LOBATTO_INT    1
#define GAUSS_DIFF           1
#define GAUSS_RADAUM_DIFF    1
#define GAUSS_RADAUP_DIFF    1
#define GAUSS_LOBATTO_DIFF   1
#define GAUSS_INTERP         1
#define GAUSS_RADAUM_INTERP  1
#define GAUSS_RADAUP_INTERP  1
#define GAUSS_LOBATTO_INTERP 1

#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <sys/time.h>

#include <vector>

#include "polylib_templated.hh"

#define USE_STL_VECTOR

/* local routines */
template<typename BaseType>
BaseType    ddot (int, BaseType* , int, BaseType* , int);

template<typename BaseType>
BaseType*   dvector (int, int);

#define START_LOG(tobj) gettimeofday(&tobj, NULL);
#define STOP_LOG(tobj) gettimeofday(&tobj, NULL); \
                        printf("Total time elapsed : %.6f sec\n" , time_diff(start, end) );

double time_diff(struct timeval x , struct timeval y)
{
  double x_ms , y_ms , diff;

  x_ms = (double)x.tv_sec*1000000 + (double)x.tv_usec;
  y_ms = (double)y.tv_sec*1000000 + (double)y.tv_usec;

  diff = (double)y_ms - (double)x_ms;
  return diff/1e6;
}


using namespace polylib ;

int main(){

  typedef double BaseType ;

  struct timeval start, end;

  int np,n,i;
  BaseType sum=0, alpha, beta;
  
#ifdef USE_STL_VECTOR
  std::vector<BaseType> z, w, p, d, dt ;
  z.resize(NPUPPER);
  w.resize(NPUPPER);
  p.resize(NPUPPER);
  d.resize(NPUPPER*NPUPPER);
  dt.resize(NPUPPER*NPUPPER);
#else
  BaseType *z,*w,*p,*d,*dt;
  z  = dvector<double>(0,NPUPPER-1);
  w  = dvector<double>(0,NPUPPER-1);
  p  = dvector<double>(0,NPUPPER-1);
  
  d  = dvector<double>(0,NPUPPER*NPUPPER-1);
  dt = dvector<double>(0,NPUPPER*NPUPPER-1);
#endif

  START_LOG(start)

#if GAUSS_INT
  /* Gauss Integration */
  printf("Begin checking Gauss integration\n");
  alpha = -0.5;
  while(alpha <= 5.0){
    beta = -0.5;
    while(beta <= 5.0){

      for(np = NPLOWER; np <= NPUPPER; ++np){
        
        GaussLegendre<BaseType> jgauss;
        jgauss.compute_points_and_weights(z, w, np, alpha, beta);

        for(n = 2; n < 2*np-1; ++n){

          jgauss.compute_polynomials(z, n, p);

          sum = ddot(np,&w[0],1,&p[0],1);
          if(fabs(sum)>EPS || sum != sum)  // if the error is greater than the tolreance or sum=NaN
            printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n" ,alpha,beta,np,n,sum);
        }
      }
      
      beta += 0.5;
    }
    printf("finished checking all beta values for alpha = %lf\n",alpha);
    alpha += 0.5;
  }
  printf("Finished checking Gauss Integration\n");
#endif

#if GAUSS_RADAUM_INT
  /* Gauss Radau Integration */
  printf("Begin checking Gauss Radau Integration\n");
  alpha = -0.5;
  while(alpha <= 5.0){
    beta = -0.5;
    while(beta <= 5.0){
      for(np = NPLOWER; np <= NPUPPER; ++np){
                
        GaussRadauLegendre<BaseType, RadauMinus> jradau;
        jradau.compute_points_and_weights(z, w, np, alpha, beta);

        for(n = 2; n < 2*np-2; ++n){

          jradau.compute_polynomials(z, n, p);

          sum = ddot(np,&w[0],1,&p[0],1);
          
          if(fabs(sum)>EPS || sum != sum)  // if the error is greater than the tolreance or sum=NaN
            printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
                      ,alpha,beta,np,n,sum);
        }
      }
      
      beta += 0.5;
    }
    printf("finished checking all beta values for alpha = %lf\n",alpha);
    alpha += 0.5;
  }
  printf("Finished checking Gauss Radau (z=-1) Integration\n");
#endif


#if GAUSS_RADAUP_INT
  /* Gauss Radau Integration */
  printf("Begin checking Gauss Radau Integration\n");
  alpha = -0.5;
  while(alpha <= 5.0){
    beta = -0.5;
    while(beta <= 5.0){
      for(np = NPLOWER; np <= NPUPPER; ++np){

        GaussRadauLegendre<BaseType, RadauPlus> jradau;
        jradau.compute_points_and_weights(z, w, np, alpha, beta);
        
        for(n = 2; n < 2*np-2; ++n){

	        jradau.compute_polynomials(z, n, p);

          sum = ddot(np,&w[0],1,&p[0],1);
          if(fabs(sum)>EPS || sum != sum)  // if the error is greater than the tolreance or sum=NaN
            printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
            ,alpha,beta,np,n,sum);
        }
      }
      
      beta += 0.5;
    }
    printf("finished checking all beta values for alpha = %lf\n",alpha);
    alpha += 0.5;
  }
  printf("Finished checking Gauss Radau (z=1) Integration\n");
#endif

#if GAUSS_LOBATTO_INT
  /* Gauss Lobatto Integration */
  printf("Begin checking Gauss Lobatto integration\n");
  alpha = -0.5;
  while(alpha <= 5.0){
    beta = -0.5;
    while(beta <= 5.0){

      for(np = NPLOWER; np <= NPUPPER; ++np){

        GaussLobattoLegendre<BaseType> jlob;
        jlob.compute_points_and_weights(z, w, np, alpha, beta);

        for(n = 2; n < 2*np-3; ++n){

          jlob.compute_polynomials(z, n, p);

          sum = ddot(np,&w[0],1,&p[0],1);
          if(fabs(sum)>EPS || sum != sum)  // if the error is greater than the tolreance or sum=NaN
            printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
            ,alpha,beta,np,n,sum);
        }
      }
        
      beta += 0.5;
    }
    printf("finished checking all beta values for alpha = %lf\n",alpha);
    alpha += 0.5;
  }
  printf("Finished checking Gauss Lobatto Integration\n");
#endif

#if GAUSS_DIFF
  printf("Begin checking differentiation through Gauss points\n");
  alpha = -0.5;
  while(alpha <= 5.0){
    beta = -0.5;
    while(beta <= 5.0){

      for(np = NPLOWER; np <= NPUPPER; ++np){

        // Jacobi<BaseType>::Utility::zwgj(&z[0],&w[0],np,alpha,beta);
        GaussLegendre<BaseType> jgauss;
        jgauss.compute_points_and_weights(z, w, np, alpha, beta);

	for(n = 2; n < np-1; ++n){
	  // Jacobi<BaseType>::Utility::Dgj(&d[0],&dt[0],&z[0],np,alpha,beta);
	  jgauss.compute_derivative_operators(d, dt, z, np) ;
	  
	  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
	  sum = 0;
	  for(i = 0; i < np; ++i) 
	    sum += fabs(ddot(np,&d[i*np],1,&p[0],1) - n*pow(z[i],n-1));
	  sum /= np;
	  if(fabs(sum)>EPS || sum != sum)  // if the error is greater than the tolreance or sum=NaN
	    printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
		   ,alpha,beta,np,n,sum);
	}
      }
      beta += 0.5;
    }
    printf("finished checking all beta values for alpha = %lf\n",alpha);
    alpha += 0.5;
  }
  printf("Finished checking Gauss Jacobi differentiation\n");
#endif

#if GAUSS_RADAUM_DIFF
  printf("Begin checking differentiation through Gauss Radau points\n");
  alpha = -0.5;
  while(alpha <= 5.0){
    beta = -0.5;
    while(beta <= 5.0){

      for(np = NPLOWER; np <= NPUPPER; ++np){

        GaussRadauLegendre<BaseType, RadauMinus> jradau;
        jradau.compute_points_and_weights(z, w, np, alpha, beta);

        // Jacobi<BaseType>::Utility::zwgrjm(&z[0],&w[0],np,alpha,beta);

	for(n = 2; n < np-1; ++n){
	  // Jacobi<BaseType>::Utility::Dgrjm(&d[0],&dt[0],&z[0],np,alpha,beta);
	  jradau.compute_derivative_operators(d, dt, z, np) ;

	  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
	  sum = 0;
	  for(i = 0; i < np; ++i) 
	    sum += fabs(ddot(np,&d[i*np],1,&p[0],1) - n*pow(z[i],n-1));
	  sum /= np;
	  if(fabs(sum)>EPS || sum != sum)  // if the error is greater than the tolreance or sum=NaN
	    printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
		   ,alpha,beta,np,n,sum);
	}
      }
      beta += 0.5;
    }
    printf("finished checking all beta values for alpha = %lf\n",alpha);
    alpha += 0.5;
  }
  printf("Finished checking Gauss Radau (z=-1) differentiation\n");
#endif

#if GAUSS_RADAUP_DIFF
  printf("Begin checking differentiation through Gauss Radau (z=1) points\n");
  alpha = -0.5;
  while(alpha <= 5.0){
    beta = -0.5;
    while(beta <= 5.0){

      for(np = NPLOWER; np <= NPUPPER; ++np){

        // Jacobi<BaseType>::Utility::zwgrjp(&z[0],&w[0],np,alpha,beta);
        GaussRadauLegendre<BaseType, RadauPlus> jradau;
        jradau.compute_points_and_weights(z, w, np, alpha, beta);

	for(n = 2; n < np-1; ++n){
	  // Jacobi<BaseType>::Utility::Dgrjp(&d[0],&dt[0],&z[0],np,alpha,beta);
	  jradau.compute_derivative_operators(d, dt, z, np) ;

	  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
	  sum = 0;
	  for(i = 0; i < np; ++i) 
	    sum += fabs(ddot(np,&d[i*np],1,&p[0],1) - n*pow(z[i],n-1));
	  sum /= np;
	  if(fabs(sum)>EPS || sum != sum)  // if the error is greater than the tolreance or sum=NaN
	    printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
		   ,alpha,beta,np,n,sum);
	}
      }
      beta += 0.5;
    }
    printf("finished checking all beta values for alpha = %lf\n",alpha);
    alpha += 0.5;
  }
  printf("Finished checking Gauss Radau (z=1) differentiation\n");
#endif

#if GAUSS_LOBATTO_DIFF
  printf("Begin checking differentiation through Gauss Lobatto points\n");
  alpha = -0.5;
  while(alpha <= 5.0){
    beta = -0.5;
    while(beta <= 5.0){

      for(np = NPLOWER; np <= NPUPPER; ++np){

        // Jacobi<BaseType>::Utility::zwglj(&z[0],&w[0],np,alpha,beta);
        GaussLobattoLegendre<BaseType> jlob;
        jlob.compute_points_and_weights(z, w, np, alpha, beta);

	for(n = 2; n < np-1; ++n){
	  // Jacobi<BaseType>::Utility::Dglj(&d[0],&dt[0],&z[0],np,alpha,beta);
	  jlob.compute_derivative_operators(d, dt, z, np) ;

	  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
	  sum = 0;
	  for(i = 0; i < np; ++i)
	    sum += fabs(ddot(np,&d[i*np],1,&p[0],1) - n*pow(z[i],n-1));
    sum /= np;
    if(fabs(sum)>EPS || sum != sum)  // if the error is greater than the tolreance or sum=NaN
      printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n",
                alpha,beta,np,n,sum);
    }
      }
      beta += 0.5;
    }
    printf("finished checking all beta values for alpha = %lf\n",alpha);
    alpha += 0.5;
  }
  printf("Finished checking Gauss Lobatto differentiation\n");
#endif

  /* check interpolation routines */
#if GAUSS_INTERP
  printf("Begin checking interpolation through Gauss points\n");
  alpha = -0.5;
  while(alpha <= 5.0){
    beta = -0.5;
    while(beta <= 5.0){

      for(np = NPLOWER; np <= NPUPPER; ++np){

        GaussLegendre<BaseType> jgauss;
        jgauss.compute_points_and_weights(z, w, np, alpha, beta);

	for(n = 2; n < np-1; ++n){
	  for(i = 0; i < np; ++i) {
	    w[i] = 2.0*i/(BaseType)(np-1)-1.0;
	    p[i] = pow(z[i],n);
	  }
	  
	  // Jacobi<BaseType>::Utility::Imgj(&d[0],&z[0],&w[0],np,np,alpha,beta);
	  jgauss.compute_interpolation_operators(w, z, d, np, np);

	  sum = 0;
	  for(i = 0; i < np; ++i) 
	    sum += fabs(ddot(np,&d[i*np],1,&p[0],1) - pow(w[i],n));
	  sum /= np;
	  if(fabs(sum)>EPS || sum != sum)  // if the error is greater than the tolreance or sum=NaN
	    printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
		   ,alpha,beta,np,n,sum);
	}
      }
      beta += 0.5;
    }
    printf("finished checking all beta values for alpha = %lf\n",alpha);
    alpha += 0.5;
  }
  printf("Finished checking Gauss Jacobi interpolation\n");
#endif

#if GAUSS_RADAUM_INTERP
  printf("Begin checking Interpolation through Gauss Radau (z=-1) points\n");
  alpha = -0.5;
  while(alpha <= 5.0){
    beta = -0.5;
    while(beta <= 5.0){

      for(np = NPLOWER; np <= NPUPPER; ++np){

        // Jacobi<BaseType>::Utility::zwgrjm(&z[0],&w[0],np,alpha,beta);
        GaussRadauLegendre<BaseType, RadauMinus> jradau;
        jradau.compute_points_and_weights(z, w, np, alpha, beta);

	for(n = 2; n < np-1; ++n){
	  for(i = 0; i < np; ++i) {
	    w[i] = 2.0*i/(BaseType)(np-1)-1.0;
	    p[i] = pow(z[i],n);
	  }

    jradau.compute_interpolation_operators(w, z, d, np, np);
	  // Jacobi<BaseType>::Utility::Imgrjm(&d[0],&z[0],&w[0],np,np,alpha,beta);

	  sum = 0;
	  for(i = 0; i < np; ++i) 
	    sum += fabs(ddot(np,&d[i*np],1,&p[0],1) - pow(w[i],n));
	  sum /= np;
	  if(fabs(sum)>EPS || sum != sum)  // if the error is greater than the tolreance or sum=NaN
	    printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
		   ,alpha,beta,np,n,sum);
	}
      }
      beta += 0.5;
    }
    printf("finished checking all beta values for alpha = %lf\n",alpha);
    alpha += 0.5;
  }
  printf("Finished checking Gauss Radua Jacobi (z=-1) interpolation\n");
#endif
#if GAUSS_RADAUP_INTERP
  printf("Begin checking Interpolation through Gauss Radau (z=1) points\n");
  alpha = -0.5;
  while(alpha <= 5.0){
    beta = -0.5;
    while(beta <= 5.0){

      for(np = NPLOWER; np <= NPUPPER; ++np){

        GaussRadauLegendre<BaseType, RadauPlus> jradau;
        jradau.compute_points_and_weights(z, w, np, alpha, beta);

	for(n = 2; n < np-1; ++n){
	  for(i = 0; i < np; ++i) {
	    w[i] = 2.0*i/(BaseType)(np-1)-1.0;
	    p[i] = pow(z[i],n);
	  }

    jradau.compute_interpolation_operators(w, z, d, np, np);

	  sum = 0;
	  for(i = 0; i < np; ++i) 
	    sum += fabs(ddot(np,&d[i*np],1,&p[0],1) - pow(w[i],n));
	  sum /= np;
	  if(fabs(sum)>EPS || sum != sum)  // if the error is greater than the tolreance or sum=NaN
	    printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
		   ,alpha,beta,np,n,sum);
	}
      }
      beta += 0.5;
    }
    printf("finished checking all beta values for alpha = %lf\n",alpha);
    alpha += 0.5;
  }
  printf("Finished checking Gauss Radau (z=1) interpolation\n");
#endif

#if GAUSS_LOBATTO_INTERP
  printf("Begin checking Interpolation through Gauss Lobatto points\n");
  alpha = -0.5;
  while(alpha <= 5.0){
    beta = -0.5;
    while(beta <= 5.0){

      for(np = NPLOWER; np <= NPUPPER; ++np){

        GaussLobattoLegendre<BaseType> jlob;
        jlob.compute_points_and_weights(z, w, np, alpha, beta);

	for(n = 2; n < np-1; ++n){
	  for(i = 0; i < np; ++i) {
	    w[i] = 2.0*i/(BaseType)(np-1)-1.0;
	    p[i] = pow(z[i],n);
	  }

    jlob.compute_interpolation_operators(w, z, d, np, np);

	  sum = 0;
	  for(i = 0; i < np; ++i) 
	    sum += fabs(ddot(np,&d[i*np],1,&p[0],1) - pow(w[i],n));
	  sum /= np;
	  if(fabs(sum)>EPS || sum != sum)  // if the error is greater than the tolreance or sum=NaN
	    printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
		   ,alpha,beta,np,n,sum);
	}
      }
      beta += 0.5;
    }
    printf("finished checking all beta values for alpha = %lf\n",alpha);
    alpha += 0.5;
  }
  printf("Finished checking Gauss Lobatto interploation\n");
#endif

#ifdef USE_STL_VECTOR
  z.clear();
  w.clear();
  p.clear();
  d.clear();
  dt.clear();
#else
  free(z); free(w); free(p); free(d); free(dt);
#endif

  STOP_LOG(end)
  
  return 0 ;
}

template <typename BaseType>
BaseType ddot (int n, BaseType *x, int incx, BaseType *y, int incy)
{
  register double sum = 0.;

  while (n--) {
    sum += (*x) * (*y);
    x   += incx;
    y   += incy;
  }
  return sum;
}


template <typename BaseType>
BaseType *dvector(int nl, int nh)
{
  BaseType *v;
  
  v = (BaseType *)malloc((unsigned) (nh-nl+1)*sizeof(BaseType));
  return v-nl;
}
